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In a one-dimensional model for a turbulent aerosol (inertial particles suspended in a random flow) 
we compute the distributions of particle- velocity gradients and the rate of caustic formation at finite 
but small Kubo numbers Ku, for arbitrary Stokes numbers St. Our results are consistent with those 
obtained earlier in the limit Ku — > and St — > oo such that Ku 2 St remains constant. We show how 
finite-time correlations and non-ergodic effects influence the inertial-particle dynamics at finite but 
small Kubo numbers. 
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I. INTRODUCTION 

Turbulent aerosols (particles suspended in turbulent flows) are ubiquitous in nature. Examples are microscopic rain 
droplets suspended in the turbulent air flow of cumulus clouds [l[ , and the motion of dust particles suspended in the 
gas surrounding a growing star (see and references cited therein). The motion of small, non- interacting particles 
suspended in a fluid is commonly modeled by the equation of motion: 

r = v . v = j(u — v) . (1) 

Here r is the particle position, v its velocity, u(r,t) is the velocity field of the flow, and 7 is the Stokes damping rate. 
Eq. {!]) assumes that the particle Reynolds number is small, that Brownian diffusion of the particles is negligible, and 
that the inertia of the displaced fluid can be neglected. 

An important dimensionlcss parameter of the problem is the ratio of the Stokes damping time 7 -1 to the correlation 
time t of the underlying flow at small length scales. This ratio is commonly referred to as the Stokes number, 
St = (71-) -1 . When the Stokes number is small, the particles are advected by the flow u(r,t). When the Stokes 
number is large, by contrast, particle inertia becomes important, allowing the particles to detach from the flow. It 
has been observed in direct numerical simulations of particles suspended in turbulent flows that particle inertia may 
give rise to large relative velocities between aerosol particles • Large relative velocities on small length scales in 
turn imply large collision rates. It is of great importance to quantitatively understand and to parameterise the St- 
dependence of collision rates of particles suspended in turbulent flows, because collision rates determine the stability 
of turbulent aerosols. 

While the dynamics of particles advected in turbulent flows is very well understood Q , inertial particle dynamics is 
more difficult to treat, and requires approximations. One possibility is to expand around the advective limit (assuming 
small Stokes numbers) @, However, this method does not capture the occurrence of singularities in the particle 
dynamics that give rise to large relative velocities at small separations P4Tl|. These singularities occur when the 
phase-space manifold (describing the dependence of particle velocity upon particle position) folds over. In fl2j these 
singularities were identified as 'caustics' analogous to light patterns on the bottom of a swimming pool on a sunny 
day. In the fold region between caustics, the velocity field at a given point in space becomes multi-valued, giving rise 
to large velocity differences between nearby particles. 

Wilkinson and Mehlig [l(| considered a limit of the problem where the flow fluctuates rapidly. The relevant 
dimensionlcss parameter is the 'Kubo number' Ku = uqt/7] where uq is the typical size of the flow velocity, and 
77 is its correlation length. The Kubo number characterises the fluctuations of u(r,t). In the limit Ku — > 0, the 
suspended particles experience the flow as a white-noise signal, and their dynamics is 'ergodic': the fluctuations of 
u(r(t),t) (and its derivatives) along a particle trajectory r(t) are indistinguishable from the fluctuations of u(ro,i) 
at a fixed point rQ. In this case the instantaneous configuration of the flow field is irrelevant to the dynamics of the 
suspended particles. In turbulent flows, the Kubo number is of order unity. Nevertheless, this approach has in the 
past yielded important insights into the dynamics of turbulent aerosols [Til, [l3|. In the limit of Ku — > 0, St — > 00 
(such that e 2 = 3Ku 2 St remains constant) it is possible to compute the fluctuations of particle- velocity gradients 
(that characterise spatial clustering of the suspended particles), the rate of caustic formation, and the distribution of 
relative velocities. 

In fbil . fl5j it was recently shown how to compute Lyapunov exponents characterising spatial clustering of inertial 
particles at finite Kubo numbers. It was found that two mechanisms for spatial clustering compete ('preferential 
concentration' [l|| and 'multiplicative amplification' [H|). 
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This raises the question of how the distribution of particle- velocity gradients and the rate of caustic formation are 
modified at finite Kubo numbers. This question is addressed in the present paper. For a one-dimensional model of 
inertial-particle dynamics we compute the distribution of particle- velocity gradients and the rate of caustic formation 
at finite Kubo numbers. 

The remainder of this paper is organised as follows. In Section [TTl the model is introduced. The distribution of 
particle- velocity gradients is computed in Section IlIIl by means of a perturbation expansion in Ku. This expansion 
describes the body of the distribution well, but does not capture its tails (related to the formation of caustics) . The 
rate of caustic formation is computed in Section IIVI using a WKB approximation valid for small Kubo numbers and 
arbitrary Stokes numbers. As in the white-noise limit, the problem can be mapped to an escape problem, but now in 
the presence of coloured noise [l?} . Finally, conclusions are summarised in Section [V] 



II. MODEL 



Understanding inertial particle dynamics at finite Stokes and Kubo numbers is a difficult problem. In this paper 
we therefore analyse inertial particle dynamics in one spatial dimension. In this case Eq. ([T]) takes the form: 



v, V = j(u(x(t),t) — v) . 



(2) 



To simplify the problem further, we take u(x,t) to be a random function with correlation time r, correlation length 
77, and typical fluctuation size Uq. We write u{x,t) = V<j){x,t) where is a homogeneous random function with zero 
mean and correlation function 



{4>{x,t)<K0,0)} =u' r) 



2 2-x 2 /(2 V 2 )-\t\/r 



(3) 



3Ku St remains 



This one-dimensional model was analysed by Wilkinson and Mehlig [T(| who computed the distribution of velocity 
gradients z = dv/dx and the rate of caustic formation in the limit Ku — s- 0, St —> 00 so that e 2 
constant. In this limit, z satisfies the Langevin equation 

i = j(A — z) — z 2 



(4) 

where A = du/dx is Gaussian white noise with zero mean and correlation function (A(t)A(t')) = 2 r ye 2 S(t — t'). The 
steady-state solution of the corresponding Fokker-Planck equation for the distribution of z is [1C| : 



P(z) 



J 

7 3 e 2 



-S(z) 



dz'e 3 ^ 



(5) 



where S(z) = z 3 /(3^ 3 e 2 ) + z 2 /(2-f 2 e 2 ) and J < is a constant probability current towards negative values of z. 
Caustics occur as z passes from —00 to +00 in a finite time. This happens at a rate determined by \ J\. The rate of 
caustic formation is given by: 



B = ±. Im 

7 2-7T 



Ai'(y) 



1 



-l/(6e 



•y=(-l/(8e 2 )) 2 /3 



/2tt 



(6) 



where the last expression is valid for small e. Eq. ^ is equivalent to the formula derived by Wilkinson and Mehlig 
fTo| . As caustics form, excursions to large values of \z\ result in slowly decaying tails of the distribution of z. It is 
easily seen that the distribution ([5]) exhibits power-law tails on the form P ~ C/\z\ 2 for large values of z. 



III. DISTRIBUTION OF z AND A AT FINITE KUBO NUMBERS 



In this section we show how to compute the steady-state distributions of particle- and flow-velocity gradients 
(z = dv/dx and A = du/dx respectively) at finite Kubo numbers. The method is based on a perturbative calculation 
of the moments of z and A. It works well when \z\ and \A\ are not too large, but fails in the tails of the distribution. 
The tails of the distribution of z are due to the formation of caustics, not described by the perturbation theory 
employed in this section. The formation of caustics at finite Kubo numbers is discussed in Section IIVI 

We introduce dimensionlcss variables: t = t'r. x = x'rj, v = v'uq, and u = u'uo, where r, 77 and uq are characteristic 
time-, space- and velocity scales discussed above. Dropping the primes to simplify the notation, Eq. becomes 



x = Ku v , i) = — (u — v) 



(7) 
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FIG. 1: Distributions of z and A. Left: difference between P(A) and Po(A). Shown are results of numerical simulations of the 
model described in Section [TT] for Ku = 0.1, St = (red o), St = 0.1 (green □), and St = 1 (blue Q). Also shown are results of 
the Ku-expansion, Eq. ()18[) extended to fourth order in Ku. Right: distribution P(z), same parameters as left. Also shown are 
the Ku-expansion Eq. ()24|) extended to second order in Ku. Also shown is the white-noise result ([5]), dashed line. 

Eq. is difficult to solve because u depends non-lincarly upon x(t). In the following we describe an approximate 
solution in terms of a perturbation expansion in powers of Ku. In [14l | this method was used to compute the Lyapunov 
exponents of inertial particles suspended in one- and two-dimensional random flows at finite Kubo numbers. 

A. Method 

The method is based on an expansion of the implicit solution of 0: 

x(t) = x + KuStoo (l - e"*/ st ) + ^ J dtt J* dt 2 e-^-^/ st u(x(t 2 ), t 2 ) . (8) 

Here xq = x(0) is the initial particle position, and vq = v(0) is the initial particle velocity. Now consider the difference 
Sx(t) = x(t) — xq between the actual trajectory of a particle and its initial position. Note that Sx(t) is proportional 
to Ku and can therefore be considered small provided Ku is small enough. In this case, one may expand u(x(t),t) in 
powers of Sx(t): 

Ou 1 d u 

u(x(t),t) = u(x ,t) + —(x ,t)Sx(t) + --^(x ,t)6x(t)Sx(t) + ... . (9) 

Inserting 5x(t) = x(t) — xq from Eq. (JSJ into Eq. (j9|) yields to second order in Ku 

dii / \ Ku Bii ( l ( ll 

u(x(t),t)=u(t)+KuStv ^(t)(l-e- t / st )+— — (t) J Ah J dt 2 e-^-^/ st u(x{t 2 )M) 

+ I K u 2 St 2 Vo 2 0(t) (l - e-^) 2 + Ku\ ^(t) (l - e-*/ st ) J* dh £ dt 2 e-^-^u(x(t 2 ),t 2 ) 

+ / dti / dt 2 dt s dt 4 e-^ +t ^-^/ st u(x(t 3 ),t 3 )u(x(U),t 4 ) + 0(Kn 3 ). (10) 

^ bt ox J Q j Jo J 
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Note that u(x(t),t) occurs on both sides of this equation. By iteratively substituting u(x(t),t) we find: 
u(x(t),t) = »(t) + KuSto„g(«) (l - e _ ' /S ') + |f jf * f' d( 2e -"'—" s '^(t) u (< 2 ) 

+ Ku2 "» I * C * e ~'"~" ),s t«>§>> (» - e - /s ') 

dti / dt 2 / di 3 / d^e-^-*^* 3 -^/ 34 — (i) — (t 2 )u(t 4 ) 



St 2 7o io Jo Jo & x ® x 

ft rt 



2 



-Ku 2 St 2 - 2 



2 ^(t) (l - e-*/ st ) 2 + Ku 2 , (l - e-*/ St ) J dfx | * dt 2 e-^)/ st J|(f)u(f 2 ) 



+ 2 S?" i, ^ J ^ Jo dt3 y o dt4e_(tl+ ' 2 " t3 "* 4)/St ^2W^ 3 ) w ^) + °( Ku3 )- ( n ) 

This equation expresses the velocity field u(x(t),t) at the particle position x(t) in terms of u(t) = u(xo,t) and its 
spatial derivatives. In the same manner higher-order contributions in Ku can be included. Eq. constitutes 
an expansion of u(x(t),t) in powers of Ku. The products of the form • • ■ d 2 u/ dx 2 (t)u(t^)u(t4) ■ ■ ■ occurring in the 
integrands on the right-hand side of Eq. can be averaged using the known statistical properties of the velocity 
field u. 

We note that in a similar manner, one may expand the strain A = du/dx (as well as higher order derivatives of u) 
along a trajectory. The result is similar to the above, but the leftmost factors §^£(t), n = 0,1, . . . in (fTT|) are replaced 
by §pr(t). Similarly particle- velocity gradients z — dv/dx can be averaged along particle trajectories. 

As outlined in the introduction, the subject of this paper are steady-state distributions of z and A. Below they are 
computed via steady-state average of moments of z and A. We denote the steady-state average of a quantity F by 

F(x(t),t)= lim i f dtF(x(t),t) = lim I f dt(F(x(t), t)) . (12) 

In the limit of T — > oo, a long trajectory can be viewed as a concatenation of many long trajectories which give rise 
to an ensemble (the average over which is denoted by (•••)) over initial configurations {x(0), v(0), u(0), <9* w(0), . . . }. 
In this limit, we expect all information about the initial configuration to be lost. At this point we therefore set 
x(0) = v(0) — in (fTTj) . We also set the initial velocity (and its spatial derivatives) to zero, it(0) = 0. That the 
solution becomes independent of the initial configuration (and hence independent of the distribution from which we 
draw the initial configuration) is explicitly shown for a number of examples in [l8j . 

In summary, the expansion described above allows to compute moments of u, v, z = dv/dx, and A = du/dx. Up to 
this point no assumption is made about the statistical properties of the velocity field u. In the following we take u to 
be Gaussian distributed; this allows us to use Wick's theorem in evaluating the terms in the perturbation expansion. 
For the particular model described in section [TTI we have (for m, n = 0, 1, 2, . . . ): 

^(x(0), tl )pM0), t 2 )\ = ( (-D ( "- m)/2 (- + n + l)l!e-l*-*l if m + n even , 
dx m dx n I [0 otherwise. v ' 

The method described above was used to compute the Lyapunov exponents for inertial particles suspended in random 
flows in one and two spatial dimensions in jl4| . In the following we show how to compute the distributions of z and 
A. 



B. Moments and distributions of z and A 

Using the method described in the previous subsection, we obtain the following expressions for the steady-state 
averages of the moments of A (for n = 0, 1, . . . ): 



A(x(t),t) 



2n+l 



jn+1 



(2n+ 1)!! 



Ku 
1 + St 



A(x(t),t) 2n = T(2n- 1)!! 



1 + 3n 



Ku 2 (l + 3St) 
(l + St) 2 (l + 2St) 



(14) 
(15) 
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The distribution P(A) = P(A(x(t), t)) is found by relating the moments of A to the Fourier transform P{a) of P(A): 



dAA n P{u,A) 



0" 



(-i)" da* 



du 



dAe 



-iaA 



P(A) 



a=0 



2ir d n 
(-i)" da' 



-P{a) 



(16) 



a=0 



Inserting this expression into the Taylor expansion of P(a) and using the moments in (|14[) and (|15[) we find 



P(a) 



x - a" 



0" 



n=0 



i\ da" 



P(a) 



a=0 



27T ^ 

n=0 



1 

2^ 



r Ku 

1 + St 



9a 2 



Ku 2 (l 



3St) 



2(l + St) 2 (l + 2St) 



-3a 2 /2 



Finally, the inverse transform yields the desired distribution: 



P(A) 



1 - 



AKu (A 2 -3)Ku 2 (l + 3St) 



1 + St 2(l + St) 2 (l + 2St) 



,-A'/o 



(17) 



(18) 



In the limit of Ku — > (or St — s- oo) this distribution converges to Pq(A) = exp[— A 2 /6]/-\/67r, the distribution of 
flow- velocity gradients evaluated at fixed position xq. At finite values of Ku, by contrast, the distribution of A is 
shifted towards negative values of A (to first order in Ku). To this order, the mean of A is 



.4 



-3Ku 
+ St 



(19) 



The fact that A is not zero is a consequence of preferential concentration [14[ . This average was computed by Wilkinson 
[TsT ] in the advective limit (St = 0) and used to determine the maximal Lyapunov exponent for inertial particles in 
one spatial dimension at finite Kubo numbers. A corresponding average figures in a one-dimensional model where the 
fluid- velocity gradients A fluctuate according to a telegraph process fl9f . 

Moreover, to second order in Ku the width of the distribution P{A) increases. The difference between P(A) and 
Po(A) is shown and compared to results of numerical simulations of the model described in Section [TT| in Fig. [I] 

We now turn to the distribution of z. Motivated by the form of (fl~8"| we try the following ansatz for P(z): a 
polynomial in z times the Ku = 0-distribution Pq(z). For small values of Ku and \z\, Po{z) is approximately Gaussian 
with zero mean and variance 3/(1 + St). From the following expressions for the steady-state moments of z 



-3Ku, 



-9Ku 



6 + 33St + 39St 2 + 10St J 



(l + St) 2 (2 + St)(l + 2St) ' 



-135Ku 



6 + 45St 



57St 2 



14St J 



l + St) 3 (2 + St)(l + 2St) ' 



-8505Ku- 



2 + 19St + 25St 2 + 6St d 



(l + St) 4 (2 + St)(l + 2St) ' 
we can determine the coefficients of the polynomial multiplying Pq(z). This gives to first order in Ku: 

2 + St-2St 2 KuSt(l + St)(6 + 9St + 2St 2 ) 3 



P(z) 



1 + St 

67T 



i (l+St)/6 



1 - Ku 



When St = 0, the distribution P(z), Eq. 
the distribution P(z) becomes {z, = KuStz) 



P(z) 



(2 + St)(l + 2St) 9 
, is identical to P{A), Eq. ((T 



(20) 
(21) 

(22) 

(23) 

(24) 



V2i 



3e 2 



-5 2 /(2 £ 2 ) 



(2 + St)(l + 2St) 
as expected. In the white-noise limit, 

(25) 



where e 2 = 3Ku 2 St. This expression is, up to a normalisation factor, identical to a series expansion of ([5]). 

Eq. (|2~4"]) is compared to results of numerical simulations in Fig. [T] The comparison shows that Eq. (f2"4"]l accurately 
describes the distribution of z at finite but small Kubo numbers provided z is not too large. The expansion employed 
here is an expansion in powers of Ku and thus assumes small fluctuations of the random variables (z, A,. . .). The 
large-z behaviour, in particular, cannot be described by this method. This implies that the algebraic tails of P (z) 
(related to the formation of caustics) arc not captured. The rate of caustic formation at finite Kubo numbers is 
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computed in the next section, using a WKB approximation. The results of the present section indicate how Eq. J6]) 
is expected to be modified at finite Kubo numbers. In the white-noise limit, where Eq. ^ applies, the rate of caustic 
formation takes an activated form [see Eq. ([6])] when e 2 = 3Ku 2 St is small. This is a consequence of the fact that the 
rate of caustic formation is determined by the rate of escape of z from z = 0, as shown by Mehlig and Wilkinson plj . 
The typical size of the fluctuations of z in the vicinity of z = is of order e 2 , and Eq. ([6]) assumes the well-known 
Arrhenius form obtained from Kramers' theory. At finite Kubo numbers, the results of this section show that the size 
of the fluctuations of z around z = are of the order of Ku 2 St 2 /(l + St). One therefore expects that the activated 
rate of caustic formation is given by (taking J to be positive) 

I ~ e~ c /( Ku2st2 ) (26) 
7 

for small values of St. In the next section we show that this expectation is borne out for very small values of St and 
compute the constant C (it turns out to be 1/96). We also compute how the small-St behaviour (f26|) rapidly crosses 
over to the large-St behaviour 

Finally, we add that in a similar fashion, the joint distribution of z and A may be obtained: 

p( z , A ) = C 1 + St) e -(A*-2A Z +(l+St)^)(l+St)/(6St) L + | _ lg StA + Q( _ 2 + 3gt + 4St 2 )z (2?) 

fi7rvSt 9(2 + St) (1 + 2St) L 



67rVSt 

- (1 + St) 2 (14 + 13St + 2St 2 )z 3 + 2A 3 (l + St) - 6(1 + St)(2 + St)^ 2 z + 3(1 + St)(8 + 9St + 2St 2 )Az 2 } 

To conclude this section we note that to first order in Kubo number, and for a Gaussian velocity field u with 
exponentially decaying time correlations (as described in section [TTJ) , the joint distribution of z and A can be obtained 
from a Langevin model: 

z = (A — 2) /St — Kuz 2 , A= -(A-A) + F . (28) 

Here F is white noise with zero mean and correlation function (F(t)F(t')) = Q5(t — if). The Ornstein-Uhlenbeck 
process for A generates the exponentially decaying time correlations. Note that A = — 3Ku/(l + St) is added to the 
fluid velocity gradients in order to keep the dominant effect of preferential concentration. We note that the joint 
distribution of z and A (|27|) satisfies the Fokker-Planck equation corresponding to Eq. (|28|) to first order in Ku. 

Eq. (|2"8")l represents the one-dimensional problem of describing the dynamics of a variable (z in this case) subject to 
coloured noise in terms of a two-dimensional white-noise problem. This approach is described in detail by Wilkinson 
[20j who solved the corresponding Fokker-Planck equation in perturbation theory. 



IV. RATE OF CAUSTIC FORMATION AT FINITE KUBO NUMBERS 



In this section we compute the rate of caustic formation at finite Kubo numbers. The formation of caustics is closely 
associated with the power-law tails of the distribution P{z). Below we employ a WKB approximation to determine 
the rate of caustic formation. For this calculation it is convenient to define dimensionlcss variables that are slightly 
different from those used in section Unl t = i/7, x = xr), v = vr)j, and u = urjj. These dimensionlcss variables were 
employed in (2l| and subsequent work on the white-noise limit of the problem. In the remainder of this section wc 
drop the tildes for ease of notation. 

Our treatment starts from (|28[1 which in the new dimensionless variables takes the form: 

i = A-z-z 2 , A = -St(A-A) + F (29) 

where F is white noise with correlation function (F(t)F(t')) = 6Ku 2 St 3 <5(i — if). A comparison between Eqs. (|28|) 
and (|29l) shows why it is appropriate to use different dimensionless variables: in this section we expand the dynamics 
in terms of small Ku in small noise levels, whereas in Sec. IIIII we expand the dynamics in terms of small Ku in the 
universal contribution due to caustics (the z 2 -term in ([28])). 

In the following we omit the non-ergodic correction A because we found that its contribution to the rate of caustic 
formation at the lowest order in Ku is negligible. In the absence of noise (F = 0), the dynamics (|2"9"]l has two fixed 
points. Their coordinates are (setting A = 0): (z^A*) = (0,0) and (z^A^) = (—1,0). The fixed point (0,0) is 
stable. In the presence of noise it becomes unstable. But for small noise amplitudes, the variables z and A fluctuate 
predominantly in the vicinity of (0,0). This corresponds to the situation described in Section IlIII However, in the 
presence of noise the variable z may escape to —00 via the unstable fixed point (—1,0). As mentioned in Section [TT1 



FIG. 2: Left: shows Hamiltonian trajectories in the z-A-plane for St = 0.2, starting at the fixed point 1. The optimal trajectory 
is shown in blue. The dashed line shows the curve A — z(z + 1). Right: same, but for St = 2. Here the dashed line shows 
A = 2z(z + l). 



this corresponds to the formation of a caustic. The rate of escape from the fixed point (0, 0) gives the rate of caustic 
formation. We note that the Langevin equation is of the form considered by Bray and McKane [l7|. Within 
a WKB approximation Bray and Kane computed the Kramers escape rate from a potential well in the presence of 
coloured noi se ( with finite correlation time). The WKB method we outline in the following is equivalent to that 
employed in [17J . 



A. WKB approximation 

The Fokker-Planck equation corresponding to (|29p is: 

r)P f) r if) f) 2 P 

^ - 5 [(, + ,» - A)P] + St-[AP] + 3KuW^ . (30) 

In the steady state we have OP/ dt = and we seek a solution of the steady-state Fokker-Planck equation of the form 

P(z, A) = exp[-S(z, A)/Ku 2 + higher orders in Ku] (31) 

(for a comprehensive description of the WKB method discussed in the following paragraphs see Dykman et al. p2j 
and references cited therein). The function S(z,A) is referred to as the 'action'. It is defined such that S = at the 
fixed point (0,0). We expect S to be a quadratic function of z and A in the vicinity of this point, corresponding to 
the case considered in Section Hill Here by contrast we are interested in the tails of the distribution, corresponding 
to large deviations of (z,A) from (0,0). Inserting the ansatz (j3"Tj) into the Fokker-Planck equation and expanding in 
powers of Ku, one obtains a first-order partial differential equation for S: 

( ^ (1 + 2)) »_ s <4f + 3S t 3(ff) J = 0. (32, 

This equation has the form of a Hamilton- Jacobi equation H(q,p) = 0, with coordinates q = (z, A) T and 'momenta' 
p = (p z ,pa) t with p z = dS/dz and pa = dS/dA. We write the 'Hamiltonian' in the standard form (see [23| and 
references cited therein): 



with 



H=p T v + \p T T>p, (33) 



- (.if +z) )• - »-(!!!«.)■ («) 



The solution of Eq. (|32[) is found by solving the Hamiltonian dynamics corresponding to (1331) : 



q = dH/dp, p=—dH/dq 



(35) 
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z* A* 


Pz 


Pa 


1 
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-1 
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1 1 
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i 




2 4 


24St 


24St 2 



TABLE I: Fixed points of the Hamiltonian dynamics (13511 . 



These equations exhibit three fixed points (q*,p*) given in Table Q] The first one (q*,Pi) = (0,0,0,0) corresponds 
to the stable fixed point of the noise-less dynamics discussed above. The second fixed point (ql,pl) = (—1,0,0,0) 
corresponds to the saddle of the noise-less dynamics. Consider solutions of (|35|) satisfying 

l a ^1 1 as t — > — oo , q(t) — > q as t — > oo , and H(q 7 p) = . (36) 
p[t) p 1 J 

To every such solution corresponds an action S^q): 

/oo 
dtpq. (37) 
-oo 

Freidlin and Wentzell [24[ gave a variational principle for the most likely path of escape satisfying the boundary 
conditions (|3f)|) . In the limit of small values of Ku, the probability distribution (|31[) is dominated by the escape path 
with extremal action. 

Fig. [2] shows numerical solutions of Eqs. ([33]) with the boundary conditions (|3"6"|) for two different values of the 
Stokes number. Shown are trajectories in the z-A-plane. These trajectories must leave the fixed point (q* 7 pl) along 
one of its unstable directions. These are found by linearising Hamilton's equations in the vicinity of (q*,Pi)'- 

'. = ( A o'-at)' M> 

Here Ai is the linearisation of the noise- less dynamics q — v with elements Aij — dvi/dqj evaluated at (<7*,p|): 

A, = ( - 1 _\ t ) . (39) 
The linearised Hamiltonian dynamics in the vicinity of (q*,p*) satisfies 

Sp = C-Hq with C- 1 = f 1 + x St - 1 \ , (40) 



and gives rise to the action 



S(q) » ^^C"^ = ^(^ 2 - 2Az + (1 + St)z 2 ) . (41) 



We see that this result is consistent with Eq. P?)) . and note that the additional factor Ku 2 St 2 is a consequence of the 
fact that different dimcnsionless variables are used in Sections Mil and IIV1 

Eq. (|4T|) [or alternatively Eq. (|40|) ] constrains the initial conditions of the Hamiltonian dynamics, leaving only 
one parameter to be varied, the angle in the z-A-planc infinitcsimally close to the first fixed point. Fig. [2] shows the 
corresponding families of trajectories. The 'optimal escape path' with extremal action is shown as a blue solid line. 
Its form (and the corresponding value of the action) depends upon the Stokes number. For large Stokes numbers, the 
optimal path approaches the curve 

A = 2z(l + z) (42) 

in the z-A-plane. In the limit of St — > oo, the z dynamics thus approaches z = z + z 2 . This is the equation of motion 
determining the optimal escape path in the white- noise limit |2l|. In the limit St — > 0, by contrast, the optimal path 
approaches the curve given by 



A = z(l + z). 



(43) 
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FIG. 3: Shows the action S in (|44|l as a function of St (solid red line). Also shown are the limiting behaviours (|50[) for large 
values of St and (|53[) for small values of St (dashed lines). Finally, results of numerical simulations of the model described 
in Section |ll] are shown (symbols). In order to allow for a quantitative comparison with the WKB theory, only parameter 
combinations that give rise to J/7 < 10 -3 are shown. Parameters: Ku = 0.01 (red,o), 0.02 (green, □), 0.03 (blue,0), 0.04 
(magenta,A), 0.05 (cyan,V), 0.07 (yellow,*), 0.1 (black,*), 0.12 (red,small o), 0.15 (green,small □), 0.17 (blue,small Q), 0.2 
(magenta, small A), 0.5 (cyan,small V), and 1 (yellow, small *). 



This condition corresponds to z = (in terms of the dimensionless units adopted in Sec. HID the condition (|43j) 
corresponds to advection, z = A). Given the optimal value S of the action, the rate of caustic formation is given by 

- ~ e" s / Ku2 . (44) 

7 

Our results for the action are summarised in Fig. [3] This figure shows results for — Ku 2 log J/7 (symbols) obtained 
from numerical simulations of the model described in Section [XXJ When J/7 is small, this expression is approximately 
given by the action. In the figure, results for different Stokes and Kubo numbers are plotted, subject to the condition 
J/7 < 10~ 3 . Also shown is the action corresponding to the optimal escape path found by numerically integrating 
(|35p. solid red line. We observe good agreement for Kubo numbers up to Ku ~ 1, despite the fact that the WKB- 
approximation is a small-Ku approximation. The asymptotic behaviours for small and large values of St (shown as 
dashed lines in Fig. [3]) are discussed in the following Subsection. 



B. Perturbation theory 

We have found an analytical expression for the action S(q) which turns out to give the correct escape action for St 
larger than St ~ 1. This expression is derived by expanding the action around q = 0: 

00 i 

S(q) = S {l \q), where = ^ af^ H , (45) 

where a)j are the expansion coefficients to order i. Provided the series (|45p converges, we may determine the action 
at the second fixed point (and thus the action of the trajectory escaping to z = —00) by evaluating (|4"5)) at z = — 1 
and ^4 = 0, i.e. S = X)So( — ^-T a i • From Eq. (|4"Tj) we know the action close to z = A = which immediately gives 
S^(q) for i < 2. We determine the higher orders S^(q) by inserting (|45|) into (|3"2"1) . requiring that all terms of order 
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i vanish. In this way we obtain the recursion 



a j+i 



J 4 

for j = 0, . . . , i — 1. Together with 



af (j + St(i - j)) + afzpU 1) - 3St 3 £ £ a£T n m+2) a<r>(i -j - m + n+ 2)(m - n) 

m— 2 n— 



(46) 



rn — 1 



a«i + a<tf\i 1) - 3St 3 ^ ^ a!r„ m+2) 4 m) (-- + n + 2)(m - n) = (47) 



this recursion determines all coefficients to order i. To third order we find 



,0) - 



fl - nQ ,3 



2(1 + St) (3) 2(1 + St) 



9St 3 (2 + St)(l + 2St) ' 1 3St 3 (l + 2St)' 

(3) 8+17St + llSt 2 + 2St 3 (3) (l + St) 2 (14+13St + 2St 2 ) 

a, = ^ , a\ = = . (48) 

2 3St 3 (2 + St)(l + 2St) 3 9St 3 (2 + St)(l + 2St) 

From these coefficients may be calculated. 

A number of comments are due at this point. First, we note that expanding the action obtained in this way in 
powers of St -1 corresponds to the small-r expansion of the escape action obtained by Bray and McKane (l7j . 

Second, expanding to lowest order in z 3 A 1 we find 

P(z, A) ~ e -( s(2,+s(3, )/ Ku2 ~ e - 5<2) / Ku2 (1 - S ( V /Ku 2 ) . (49) 



The last expression corresponds to the small-Ku limit of (|27|) . Factors of KuSt arises because different dimcnsionlcss 
units are used in Section IIIII and here. 

Third, expanding the action to + at z = — 1 and A = to lowest order in St -1 results in S = 1/(18 St). 
This is again the white-noise result ©• 

(i) — l 

Fourth, expanding the higher-order coefficients in St gives 

1 1 1 4 
S = — — + 5 r + 7 + . . . . 50) 

18 St 90 St 3 105 St 5 189 St 7 V ; 

This result can also be obtained from the lowest order of the small-r expansion by Bray and McKane [13] (Eq. (12) 
in their paper). The lowest order of this expression is the white- noise result, shown as a dashed line in Fig. [3] 

The small-St asymptotics of the action can be derived following a procedure outlined in (l7j . The starting point 
is Eq. (|43[) . This condition can be interpreted as a fixed-point condition for the z-dynamics which rapidly adjusts to 
the slowly- moving variable A. This fixed point becomes unstable for z < — 1 (and for z > 0). In the former case, 
the z-coordinate escapes from z = —1 to z = — oo with A = 0. Solving Eq. (|43l) for z in the range —1 < z < gives 
z^ = — 1/2 ± \/l + 4A/2. The solution z^' is appropriate in the interval — 1 < z < —1/2, while the solution 
pertains to —1/2 < z < 0. There are two corresponding solutions of Hamilton's equations subject to the condition 

(USD 

pW = , = 4 +) e St 7(3St 2 ) , AM = A«e s « , (51) 

p(-)=0, pM=0, A(-)=4- ) e- stt . (52) 

The first solution (denoted (+)) describes how A decreases from a small negative value Aq + ^ at z = to A = —1/4 
at z = —1/2. Then the second solution (denoted by (— )) takes over, describing how A increases to A = at z = — 1 
asM oo. The action vanishes for the (— )-branch, and the action for the (+)-branch gives: 

S = — 7t . (53) 
96 St 2 

This result is equivalent to the lowest order in a large-r expansion given in [l7| . Eq. (|53[) is shown as a dashed line 
in Fig. [3j Our numerical results (solid red line) converge to this asymptote, albeit very slowly. 
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V. CONCLUSIONS 



In this paper we have computed the distribution of particle- velocity gradients z = dv/dx and the rate of caustic 
formation at finite Kubo numbers, for a one-dimensional model of a turbulent aerosol. The body of the distribution 
was computed by means of a perturbation expansion in Kubo number, and the rate of caustic formation within a 
WKB approximation. We have described the differences between these two approaches, and how they arc consistent. 
The results complement and extend results obtained earlier in the white- noise limit Ku —> and St — > oo, such that 
e 2 = 3Ku St remains constant. These earlier results have yielded important insight into the mechanisms governing 
inertial particle dynamics, but it is of interest to ask how these results are modified at finite Kubo numbers - since 
turbulent flows have Kubo numbers of order unity. 

A finite Kubo number has two effects: the fluid velocities have a finite correlation time, and the suspended particles 
may explore configuration space preferentially. Both affect the distribution of particle-velocity gradients and the rate 
of caustic formation. To lowest order in Ku, modifications are due to the finite correlation time of the flow. The effect 
due to preferential sampling occurs at higher order in Ku. 

We find that the rate of caustic formation exhibits an activated St-dependence, as in the white-noise limit. For 
St > 1, our result for the rate of caustic formation agrees with earlier results obtained in the white- noise limit: 

I „ c -l/(6 £ 2 ) with £ 2 = 3Ku 2 St _ ( 54 ) 

7 

For very small Stokes numbers (St < 0.04) we find, by contrast, that 

- ~ e-V0*Ku»st») . ( 55) 
7 

This scaling is consistent with a parameterisation of the collision rate for particles in turbulent aerosols proposed in 
0. At intermediate Stokes numbers (0.04 < St < 1 for the model considered here) the St-dependence of the action 
is more complicated, but the WKB approximation yields a quantitative description. The action for the caustic-rate 
formation is shown in Fig. [3l The crossover between the two asymptotic expressions (|54| and ([55]) occurs at St = 3/16 
for the model considered here. 

While we expect essentially similar results in two and three spatial dimensions it is nevertheless important to check 
this in detail. A difference between one and higher dimensions is of course that one-dimensional flows are always 
compressible. Using the approach described in Section |n] it will be possible to compute the probability that particles 
are found in vortical regions in two- and three-dimensional incompressible flows at finite Kubo numbers. This is of 
interest since this probability was measured in direct numerical simulations of inertial particles suspended in turbulent 

flows m. 

In one spatial dimension, a number of open questions still remain. First, the WKB approximation was only 
performed to lowest order in Ku, the prefactor of the exponential in the distribution of z must also be computed. It 
is expected that it will give rise to power-law tails ~ \z\~ 2 of the distribution of z. Second, it remains to be seen to 
which extent preferential sampling at finite Kubo numbers modifies this prefactor. 

Acknowledgements. We gratefully acknowledge financial support by Vetenskapsradet and by the Goran Gustafsson 
Foundation for Research in Natural Sciences and Medicine. 
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